HABITAT SELECTION MODELS FOR PACIFIC SAND LANCE (AMMODYTES HEXAPTERUS) IN PRINCE WILLIAM SOUND, ALASKA

Abstract We modeled habitat selection by Pacific sand lance (Ammodytes hexapterus) by examining their distribution in relation to water depth, distance to shore, bottom slope, bottom type, distance from sand bottom, and shoreline type. Through both logistic regression and classification tree models, we compared the characteristics of 29 known sand lance locations to 58 randomly selected sites. The best models indicated a strong selection of shallow water by sand lance, with weaker association between sand lance distribution and beach shorelines, sand bottoms, distance to shore, bottom slope, and distance to the nearest sand bottom. We applied an information-theoretic approach to the interpretation of the logistic regression analysis and determined importance values of 0.99, 0.54, 0.52, 0.44, 0.39, and 0.25 for depth, beach shorelines, sand bottom, distance to shore, gradual bottom slope, and distance to the nearest sand bottom, respectively. The classification tree model indicated that sand lance selected shallow-water habitats and remained near sand bottoms when located in habitats with depths between 40 and 60 m. All sand lance locations were at depths <60 m and 93% occurred at depths <40 m. Probable reasons for the modeled relationships between the distribution of sand lance and the independent variables are discussed.

Pacific sand lance (Ammodytes hexapterus) play an important ecological role as energyrich prey (Anthony and Roby 1997) for seabirds, marine mammals, and predatory fishes in Prince William Sound (PWS), Alaska (Kuletz and others 1997). Elsewhere, sand lance population dynamics have been correlated to the re-86(3) row in sandy substrates, avoiding rocky, muddy, and coarse gravel bottoms (Reay 1970). By burrowing into substrates for refuge (while not foraging or during winter) sand lance effectively reduce their risk of predation and energy expenditure (Pinto and others 1984). Sand lance burrow in both subtidal (Hobson 1986) and intertidal (Dick and Warner 1982) substrates and are associated with moderately sloped beaches composed of well-washed gravel or finer sands (Dick and Warner 1982). Wright and others (2000) also observed that lesser sand eels (A. marinus) associate with sandy substrates and avoid those with Ͼ10% fine material. While not burrowed, Pacific sand lance generally remain close to sandy refuges (Kò hlmann and Karst 1967;Hobson 1986). Kò hlmann and Karst (1967) observed that sand eels (Ammodytes sp.) of the Baltic Sea foraged diurnally, moving farther offshore during daylight hours and returning to burrow in shallow nearshore sand bottoms at night.
We infer from the apparent ecological importance of sand lance that resource managers and planners will find it useful to be able to predict important sand lance habits in order to maintain productivity of nearshore ecosystems. The studies above present descriptive information on habitat variables associated with sand lances (Ammodytes spp.). We have accepted that each of the variables presented above are important factors in sand lance habitat selection. We did not seek to find new habitat associations; rather, our goals were to gain further understanding of the relative importance of these variables and to demonstrate how they might be used to model Pacific sand lance habitat selection.

Study Area
We conducted this study in PWS, an embayment of about 10,000 km 2 located on the southcentral coast of Alaska. The climate is maritime with a mean annual precipitation of approximately 1.6 m and moderate air temperatures for the subarctic. The coastline of PWS is rugged, with mountains Յ4000 m in elevation, numerous islands, fjords, and tidewater glaciers. Nearshore bathymetry is characterized by both shallow water shelves and steeply sloping bottoms.
Four study areas were selected in PWS: (1) the northern study area, which included Valdez Arm and Port Valdez; (2) the central study area, which included waters near Naked and Knight Islands; (3) the southern study area, which included Icy and Jackpot Bays; and (4) the north shore of Montague Island (Fig. 1). These locations were selected because seabird foraging studies indicated that sand lance were present within each area (Irons and others 1997;Suryan and others 2000).

Hydroacoustic Data Collection
Some habitat variables (depth, bottom slope, and bottom type) were determined through the analysis of hydroacoustic data. Haldorson and others (1997) conducted both nearshore and pelagic hydroacoustic and catch studies within PWS. Sand lance were only detected nearshore, consistent with the observations of Raey (1970). Therefore, we chose to utilize the most extensive nearshore hydroacoustic data set, which was collected 17-27 July 1997 (Haldorson and others 1997;Ostrand and others 2005). This survey used a commercial purse seiner, the 18m F/V Miss Kaylee, as a data-collection platform. Data were collected with a single beam 120 kHz BioSonics DT4000 system with a 6Њ beam angle. Transects were run at 11 kph with the transducer towed beside the vessel. The effective range of the equipment was 117 m from the transducer. Location data were obtained from a precision lightweight global positioning receiver (PLGR). PLGR units have a worst-case horizontal position accuracy of 10 m at speeds Ͻ36 kph (Anonymous 1995).

Hydroacoustic Survey Design
The hydroacoustic study used 1-ϫ 12-km shoreline blocks located within the 4 study areas (Haldorson and others 1997;Ostrand and others 2005). Prior to the survey, contiguous blocks were delineated in study areas, which included all available shoreline. Due to the large extent of the northern and southern areas and the impracticality of sampling their entire shorelines, alternate blocks and 1 additional randomly selected block were deleted from these areas. All possible blocks were retained in the central and Montague areas. The ultimate study design contained 9, 8, 8, and 2 blocks in the north, central, south, and Montague Island areas, respectively (see Table 1  data on extent of survey). Within each block, 20 continuous 1.2 km transects were laid out in a zigzag pattern (Fig. 2). The nearshore apexes of transects were located as close to shore as possible.
Semivariograms produced by a spherical kriging algorithm (Oliver and Webster 1990) suggested that hydroacoustic transect data were spatially auto-correlated to a distance of 0.5 km. We then delineated a 0.5-km buffer on both sides of the hydroacoustic transects to define the extent of habitat to be sampled and modeled (Table 1; Fig. 2).

Sand Lance Locations
A set of 29 sand lance locations was derived from collections and observations that were recorded by studies in PWS during the summers of 1996 to 1999 (Table 2). Identifications were made to species but not to age class. Position coordinates for all locations were determined by Global Positioning System (GPS) receivers (100 m; Leick 1992; Anonymous 1995) during the original studies. There may have been inconsistent biases associated with the sampling methods used to obtain sand lance locations. To investigate possible bias in the detection methods, we grouped the sites into 5 classes: purse seine locations obtained during hydroacoustic surveys (n ϭ 7), observations made while radio tracking black-legged kittiwakes (n ϭ 13), samples dip-netted beneath foraging marbled murrelets (n ϭ 5), SCUBA diver observations (n ϭ 3), and opportunistic observation (n ϭ 1). We tested for bias among methods by conducting ANOVAs (SAS 1996) for each continuous variable and likelihood ratio Chisquare (SYSTAT 1997) for each binary variable. The single opportunistic observation was excluded from testing for bias. For the ANOVAs and Chi-square test, we considered P Յ 0.05 significant.
We selected 58 random locations to represent habitats available to sand lance for selection. Conceptually, the random data set includes available locations in proportion to their presence within the study areas. Available locations include habitats that range from suitable to unsuitable for selection by sand lance (Manly and others 1993). We have assumed that the sand lance location data set contains only suitable habitat that has been selected by sand lance from all available habitat within the study area. Both sand lance and available locations were within the buffered hydroacoustic survey route. For all locations, we determined values for bottom type (sand or not sand), depth, shoreline type (beach or not beach), distance to shore (at mean tide), bottom slope, and distance to the nearest sand bottom (located within the buffered sampling area).

Bottom Type
To classify the seabed, we analyzed the hydroacoustic data (Haldorson and others 1997;Ostrand and others 2005) with bottom typing software (VBT Seabed Classifier, BioSonics, Inc., Seattle, WA). This software produced 6 variables that described the characteristics of the bottom signal: energy of the sediments echo (E0), energy of the 2nd part of the 1st bottom echo (E1), energy of the 2nd bottom echo (E2), energy of the 1st part of the 1st bottom echo (E1'), sediment thickness (ST), and fractal dimension of the 1st bottom echo (FD). The data set did not consistently have a 2nd bottom echo; therefore, we omitted E2 from further analysis. To examine possible interactions between parameters, we created additional variables by multiplying values for E0 and E1, E0 and FD, and E1 and E1'. We adjusted the software to average all variables over 30-m horizontal intervals.
To calibrate the bottom typing, during summer 1998, we obtained sediment samples with a Ponar grab at randomly selected sites located on the hydroacoustic track at depths Ͻ100 m.
The sampling device was dropped from a boat within Ϯ10 m of the pre-selected location. Positions were determined with a PLGR. Samples Ͼ50 g, dry weight, were obtained at 20 of 49 sites at which grabs were attempted. Replicate samples were collected at 5 sites. Our rate of FIGURE 2. Paired maps of the 4 study areas depicting Pacific sand lance (Ammodytes hexapterus) and random locations, hydroacoustic transects with 0.5-km buffer, sand bottoms, beach and non-beach shorelines, and water depth.
sampling success indicated to us that the Ponar grab was not of sufficient size to consistently sample the bottom within PWS. To improve sampling success rate we changed gear the following year. During the summer of 1999, we collected samples from 28 of 30 sites with a Shipek bottom grab. Duplicate samples were collected at 23 of these sites. We also collected 21 samples from 24 sites with a bottom dredge with a 28 ϫ 76-cm opening and 76-cm depth, during July 1999. One sample was collected at each dredge site. The final bottom typing validation data set contained 69 samples (Fig. 1).
Sediment samples were frozen and then oven dried (150Њ C for 3 h) prior to laboratory analysis. Grain size analysis was performed on sediment samples using a sieve procedure (Day 1965), which determined percentage gravel, sand, silt, and clay for each sample following the USDA scale (Gee and Bauder 1986). For sites where 2 samples were taken, they were analyzed separately and the results were averaged. Samples were designated as sand, Ն80% sand, or not sand (Folk 1980).

86(3)
To ascribe the bottom types throughout our study areas we used S-plus (MathSoft 1999) classification tree modeling (Chambers and Hastie 1992;Bell 1996). We used the 69 sediment classifications of the validation data as the dependent variable and the associated values of the characteristics of the hydroacoustic bottom signals as independent variables for the training data set. Trees were overfit to the data, setting S-plus options of minimum deviance ϭ 0 and minimum size ϭ 2. Model selection and misclassification rate estimation were based on a jackknife procedure. Initially we removed the 1st observation from the data set and fit the tree to the remaining observations. The number of terminal nodes (size) of the tree was reduced (pruned) to each size between 1 and the unpruned tree size, and the 1st observation was classified by each of these trees. The pruning function uses cost-complexity to determine the optimal tree of a specified size (Chambers and Hastie 1992;Bell 1996). We noted if the classifications of the 1st observation were correct. Next, we replaced the 1st observation and removed the 2nd observation and repeated the pruning and classification procedures. This method was continued until each observation had been removed from the data set and reclassified by trees of each size. The %-misclassification rate (the number of incorrect classifications divided by the number of observations ϫ 100) was calculated for each tree. We selected the tree size having the lowest jackknifed misclassification error rate to predict the bottom type of the sampling areas. We considered misclassification rates Ͼ25% to be unacceptable.

Depth, Bottom Slope, Distance to Shore, and Distance to the Nearest Sand
The bathymetry within the sampling area was determined by kriging (Oliver and Webster 1990) depth values extracted from the hydroacoustic data set. The maximum depth recorded by hydroacoustics was 117 m. All locations within the sampling area deeper than 117 m were excluded from analysis. Sand lance are most commonly found in shallow water (Ͻ50 m) and rarely occur at depths of Ͼ100 m (Robards and Piatt 1999) justifying this demarca-tion. To account for the slope from shore, adjacent coastline depths were made equal to 0 and incorporated in the kriging. A slope coverage, expressed in degrees (0 to 90Њ), was calculated within ARC/INFO from depth and location data coverages. Distance to shore and distance to the nearest sand bottom values were determined for sand lance and random locations using ARC/INFO.

Shoreline Type
Shoreline habitat type was determined through aerial visual observations conducted during spring 1999. These data were digitized and converted into a GIS coverage (Research Planning Inc., PO Box 328, 1121 Park Street, Columbia, SC, unpubl. data). Shorelines were classified as exposed rocky shores, exposed solid man-made structures, exposed wave-cut platforms in bedrock, fine-to medium-grained sand beaches, mixed sand and gravel beaches, gravel beaches, riprap, exposed tidal flats, sheltered rocky shores, sheltered riprap, vegetated steeply-sloping bluffs, sheltered tidal flats, and salt-and brackish-water marshes. We conducted a preliminary analysis with the 2 modeling methods (see modeling techniques below) and found no relation between shoreline class and distribution of sand lance. The beach shoreline types (fine-to medium-grained sand beaches, mixed sand and gravel beaches, and gravel beaches) were descriptively similar to the intertidal habitats that Dick and Warner (1982) and Robards and Piatt (MD Robards and JF Piatt, Biological Research Division, US Geological Survey, 1011 E. Tudor Road, Anchorage, AK, unpubl. data) observed to be sand lance habitat in areas near PWS. Therefore, we collapsed shoreline types into beaches and nonbeaches for final analyses.

Habitat Selection Models
We applied both logistic regression (Manly and others 1993;SAS 1996) and classification tree (Chambers and Hastie 1992; Bell 1996; MathSoft 1999) methods to model sand lance habitat selection. Because these methods examine data differently and may produce different results (Dettmers and others 2002), we chose to use both methods to obtain a broad prospective on habitat selection. For both methods, we compared the characteristics of known sand lance locations to available (random) lo- cations. Presence of sand lance or available location was the dependent variable; bottom type, depth, distance to shore, bottom slope, distance to the nearest sand, and shoreline type were the independent variables. Presence of sand lance or available location, bottom type, and shoreline type data were binary and all other variables were continuous. Prior to conducting regressions, variables were checked for independence through correlation analysis (SAS 1996). We considered r Ͻ 0.50 to indicate independence. We then developed a model set composed of all possible combinations of 6 variables excluding interactions and higher order terms (63 models). Logistic regression was fitted to all equations within the model set, and these were ranked based upon Akaike's information criterion corrected for small sample sizes (AICc) (Akaike 1973; Burnham and Anderson 1998). The AICc statistic describes the fit of a model while penalizing for variables that explain minimal error. Burnham and Anderson (1998) have suggested that all models whose AICc values differ from the lowest observed AICc by Ͻ2 (difference from lowest AICc, hereafter ⌬ AICc) should receive consideration. We examined all models that had ⌬ AICc Ͻ2 and made descriptive comparisons of them. We also determined importance values (Burnham and Anderson 1998) from the set of all possible models for each independent variable using the candidate model set. Importance values are the sum of Akaike weights of the models in which the variable is included. Classification tree model selection followed the same procedures as presented under bottom typing.

RESULTS
The logistic regression and classification tree models all indicated a strong association between sand lance distribution and depth, with weaker and inconsistent relationships to beach shorelines, bottom type, distance to shore, bottom slope, and distance to the nearest sand bottom (Table 3, Fig. 3). Fifteen logistic regression models had ⌬ AICc Ͻ2 and all of these were highly significant (P ϭ 0.0001; Table 3). All of these models contained depth as a variable, with an individual variable P Ͻ 0.05. No other variable appeared in all 15 models. Distance to the nearest sand bottom was the only variable that did not appear in any of the 15 models. The importance values of the independent variables, out of a maximum value of 1.00, were depth ϭ 0.99, shoreline type ϭ 0.54, bottom type ϭ 0.52, distance to shore ϭ 0.44, slope of the bottom ϭ 0.39, and distance to the nearest sand bottom ϭ 0.25. The classification tree model had a jackknife misclassification of 10%.  The tree indicated that at depth Ͼ40 m, habitat selection was based upon distance-to-sand and depth. We further examined the distance-tosand variable at depths Ͼ40 m and determined that x ϭ 37.7, s ϭ 53.3 m, n ϭ 2 for sand lance locations, and x ϭ 610.9, s ϭ 741.5 m, n ϭ 36 for available sites.
Values (x Ϯ s) for continuous variables at locations where sand lance were present and at available locations were 21.1 Ϯ 16.7 and 53.1 Ϯ 36.3 m depth; 418.3 Ϯ 377.3 and 688.4 Ϯ 584.2 m distance from shore; 3.6 Ϯ 2.9 and 5.7 Ϯ 6.2 degrees slope; 528.0 Ϯ 567.6 and 590 Ϯ 735.3 m distance from the nearest sand bottom, respectively (see Fig. 4 for comparative histograms). Sand bottoms were associated with 20.7% (6) of the sand lance locations and 8.6% (5) of the available locations. Beach was the nearest shoreline type at 90.0% (26) and 69.0% (40) of the sand lance and available locations, respectively. We further examined our data to determine if there was a relationship between beach (a binary variable) and bottom slope. At available locations, bottom slopes associated with beach shorelines (x ϭ 4.2, s ϭ 4.8, n ϭ 40) were less steep than slopes associated with nonbeach shorelines (x ϭ 9.2, s ϭ 7.7, n ϭ 18; Kruskal Wallis 1-way ANOVA, H ϭ 5.28, P ϭ 0.022).
All comparisons of the methods of collecting sand lance locations yielded non-significant results (ANOVA test for depth P ϭ 0.81, slope P ϭ 0.77, distance to the nearest sand bottom P ϭ 0.60, and distance to shore P ϭ 0.09; likelihood ratio Chi-square tests for sand or notsand bottom P ϭ 0.20 and for beach or notbeach shoreline P ϭ 0.67). Correlation analysis indicated that all possible pairs of variables had r Ͻ 0.50, thereby meeting our requirements for independence for all modeling components.
The bottom typing classification tree had seven 7 nodes and utilized the following variables: energy of the 1st part of the 1st bottom echo (E1'), energy of the 2nd part of the 1st bottom echo (E1), fractal dimension (FD), energy of the sediments echo (E0) ϫ FD, and E1 ϫ E1'.
The jackknife estimated misclassification rate was 6%. We also attempted to model more than 2 bottom types; however, these attempts resulted in unacceptable misclassification rates.

DISCUSSION
The models presented indicated that sand lance distribution was strongly associated with 86 (3) shallow water depths. Few of the sand lance locations were at depths Ͼ40 m and none at depths Ͼ60 m (Fig. 4), agreeing with previous estimates of sand lance depth distribution . Selection for shallow water is consistent with the findings of Winslade (1974) that sand lance are visual foragers and are sensitive to light while burrowed. Hence, the limit of light penetration through the water column may be the mechanism behind the selection of shallow habitats.
The importance values for the substrate variables shoreline type and bottom type ranked 2nd and 3rd. These model outcomes indicate selection for beach shorelines over non-beaches and sand bottoms over other bottom types. The linkage between sand lance and sand substrates as burrowing habitat and refuge is well documented (Dick and Warner 1982;Hobson 1986), and we had expected a strong association between sand lance distribution and these parameters.
The importance values indicated that the association between sand lance distribution and distance to shore and bottom slope ranked below depth and the substrate variables. A relationship between distribution and distance to shore intuitively follows from the strong selection for shallow habitat. Shallower water is expected closer to shore; however, our correlation analysis indicated that there was not a correspondence between depth and distance to shore. Within PWS, at the spatial scale of this study, shallow water was available both close to and distant from shore, which may explain the importance ranking of the variables. Our examination of bottom slope relative to depth suggests that moderate bottom slopes are associated with beach habitats in PWS. Therefore, the inclusion of bottom slopes within logistic regression models having ⌬ AICc values Ͻ2 may be the result of covariance with beach shorelines.
The classification tree analyses yielded a model with lower misclassification rate for bottom type than for habitat selection. This may have been a partial result of the nature of the data sets examined in each analysis. Bottom typing compared characteristics of known sand locations to non-sand sites; whereas, the habitat selection model compared known sand lance locations to available, or random, sites. The presence of suitable sites within the avail-able data set may have resulted in overfitting (Bell 1996;Burnham and Anderson 1998) as the tree modeling process attempted to classify the suitable available sites as non-sand-lance locations. In examining the habitat selection tree (Fig. 3), the data are 1st split at 40-m depth. At depths Ͻ40 m are 5 terminal nodes, 1 of which conferred a non-sand-lance habitat type. This node indicated that if depth was 13 to 33 m and the distance to the nearest sand bottom was 69 to 560 m then the habitat was unsuitable for sand lance. To us, this branching was counterintuitive and appeared to be the result of the tree accommodating available locations that were suitable for sand lance, indicating overfitting. The presented tree was determined through a rigorous and objective jackknife procedure; however, in light of the above arguments we suggest a modified model that indicates that habitat selection is based solely upon depth when depth is Ͻ40 m. The resulting change would be represented by a reduction to 1 terminal node on the left side of the tree, below the initial split of the data at 40-m depth.
The habitat selection tree indicated that at depth Ͼ40 m sand lance did not venture far from sand bottoms, consistent with the observations of Hobson (1986). This relationship was not reported by logistic regression. The recursive partitioning of the data by the classification model (Chambers and Hastie 1992;Bell 1996) has revealed a pattern in the data that would have otherwise remained undetected and illustrates the value of multiple statistical approaches. However, we must caution against the over interpretation of this result due to the small number (n ϭ 2) of sand lance locations observed at depth Ͼ40 m. We suggest that this is an interesting result worthy of further study. The comparison of methods used to determine sand lance locations failed to detect biases in sampling methods. This finding substantiated our application of these data to develop the presented models. This was not a definitive comparison of data collection methods, and we caution against making further inferences from these results. The total sample size of verified sand lance locations was less than desirable to make inference to an area as large as PWS. However, these data were collected by 5 different studies over 4 y. Each of the studies had an interest in locating sand lance and put an emphasis on recording their locations. The low sample size in light of such an effort suggests that locating sand lance is a challenging task. It is improbable that a greater effort to locate sand lance in PWS will be initiated in the near future.
There are apparent biases in our study toward near shore and shallow water habitats. Our initial pilot work (Haldorson and others 1997) determined that sand lance were only located near shore, and we designed our hydroacoustic survey accordingly. Had we conducted our surveys throughout the PWS embayment, our models would have indicated a strong preference for near-shore habitats. Interpretation of models must be cognizant of the assumption of near-shore habitat selection. Assumptions that our survey was biased toward shallow water habitats are reflective of a lack of knowledge of the bathymetry of PWS. Nearshore habitats of PWS are characterized by shallow sloping beaches and precipitous cliffs that extend below the water's surface. Our correlation analysis, which found no relationship between water depth and distance from shore, is reflective of the extremes in nearshore bathymetry. Our methods of detecting sand lance are also independent of depth with the exception of SCU-BA diving, which is limited by light penetration of the water column. However, our ANOVA analysis of detection methods did not find a significant difference of depths by method, suggesting that the effective depth of SCUBA was sufficient.
Our models were consistent with earlier descriptive studies and have provided insight into the association between sand lance and pelagic habitats. We have shown that of the parameters investigated, depth is dominant followed by shoreline type and bottom type in influencing habitat selection. However, we wish to caution against inference that sand lance select habitat in the same manner throughout their range. We speculate that if the availability of habitat differs greatly from what we observed, selection may change (Manly and others 1993). Furthermore, sand lance may have a degree of plasticity in their habitat preference (for example, Pinto and others 1984). Habitat choice may be the result of complex interactions among all of these variables.
We hope that the models presented will be useful to natural resource managers and other researchers in identifying potential pelagic sand lance habitats. Each model can be used to develop GIS coverages that will indicate probable sand lance habitat, provided data on the independent variables are available. We have not attempted to select a preferred model or modeling method. We suggest that there is utility in examining multiple models and we leave to readers the choice of model to be selected for their purposes. Sand lance habitats are likely to be important to the many animals that prey upon sand lance (Willson and others 1999) and may require protection from anthropogenic disturbance and excessive development. We recommend testing the models against independent data or conducting field verification prior to accepting the predictions.